US COVID Data

Data Prep

us <- read.csv("https://raw.githubusercontent.com/JaclynCoate/6373_Time_Series/master/TermProject/Data/USdaily.7.26.20.csv", header = T, strip.white = T)
us <- transform(us, date = as.Date(as.character(date), "%Y%m%d"))
us <- subset(us, select = -c(states, dateChecked, hospitalized, lastModified, total, posNeg, totalTestResultsIncrease, hash))
us[is.na(us)] <- 0
us = us[order(as.Date(us$date, format = "%Y%m%d")),]
head(us)
##           date positive negative pending hospitalizedCurrently
## 187 2020-01-22        2        0       0                     0
## 186 2020-01-23        2        0       0                     0
## 185 2020-01-24        2        0       0                     0
## 184 2020-01-25        2        0       0                     0
## 183 2020-01-26        2        0       0                     0
## 182 2020-01-27        2        0       0                     0
##     hospitalizedCumulative inIcuCurrently inIcuCumulative
## 187                      0              0               0
## 186                      0              0               0
## 185                      0              0               0
## 184                      0              0               0
## 183                      0              0               0
## 182                      0              0               0
##     onVentilatorCurrently onVentilatorCumulative recovered death
## 187                     0                      0         0     0
## 186                     0                      0         0     0
## 185                     0                      0         0     0
## 184                     0                      0         0     0
## 183                     0                      0         0     0
## 182                     0                      0         0     0
##     totalTestResults deathIncrease hospitalizedIncrease negativeIncrease
## 187                2             0                    0                0
## 186                2             0                    0                0
## 185                2             0                    0                0
## 184                2             0                    0                0
## 183                2             0                    0                0
## 182                2             0                    0                0
##     positiveIncrease
## 187                0
## 186                0
## 185                0
## 184                0
## 183                0
## 182                0

Stationarity: Current Hospitalizations

It is difficult to assume stationarity for this data due to multiple factors. We are working under the assumption that COVID is a novel virus and cases as well as hospitalizations will eventually return to zero. This being said our current modeling techniques do things such as return to the median or mimic the previously seen trends. Also, we see a heavy wandering trend in both new cases and hospitalization would be dependent on this as well as time. We will review the data and see what, if any, non-stationary components reveal themselves and model the data accordingly.

Current Hospitalization Realization

Traits:

  • Heavy wandering behavior
  • What appears to be some noise that could be pseudo-cyclic behavior hidden by the large numbers.
ggplot(data = us, aes(x=date, y=hospitalizedCurrently))+
  geom_line(color="orange")+
  labs(title = "Current COVID Hospitalized Cases US", y = "Thousands", x = "") +
    theme_fivethirtyeight()

Independent Variable Review

When reviewing the variables and the scatter plot from our previous EDA we can see that there are some correlations that were expected, for example currentHospitalizations is correlated with the variables that refelct ICU and Ventilaor patients. These metrics wouldn’t exist without hospitalizations. These variables also cannot help up predict hospitalizations because they after occurances. So we will actually be leveraing variables such as positive increase in order to see if there is some sort of correlation between hospitalized patients and the number of positive cases.

  • Positive Increase Trend
ggplot(data = us, aes(x=date, y=positiveIncrease))+
  geom_line(color="orange")+
  labs(title = "Increase in COVID Cases US", y = "Thousands", x = "") +
    theme_fivethirtyeight()

Multivariate MLR w/ Correlated Errors for Currently Hospitalized Patients

Forecast Independent Variables: Increase Positive Cases

  1. Evaluation: Slowly dampening ACF and heavy wandering consistent with a (1-B) factor. As well as frequency peaks at 0 and .14. Consistent with (1-B) and seasonailty component of 7.
#a
plotts.sample.wge(us$positiveIncrease)

## $autplt
##  [1] 1.0000000 0.9740745 0.9439116 0.9109437 0.8869924 0.8711719 0.8605261
##  [8] 0.8462935 0.8144161 0.7772852 0.7356766 0.7045767 0.6829609 0.6646632
## [15] 0.6432819 0.6088202 0.5688872 0.5293417 0.5001353 0.4761589 0.4592724
## [22] 0.4426271 0.4148783 0.3789333 0.3458221 0.3182895
## 
## $freq
##  [1] 0.005347594 0.010695187 0.016042781 0.021390374 0.026737968
##  [6] 0.032085561 0.037433155 0.042780749 0.048128342 0.053475936
## [11] 0.058823529 0.064171123 0.069518717 0.074866310 0.080213904
## [16] 0.085561497 0.090909091 0.096256684 0.101604278 0.106951872
## [21] 0.112299465 0.117647059 0.122994652 0.128342246 0.133689840
## [26] 0.139037433 0.144385027 0.149732620 0.155080214 0.160427807
## [31] 0.165775401 0.171122995 0.176470588 0.181818182 0.187165775
## [36] 0.192513369 0.197860963 0.203208556 0.208556150 0.213903743
## [41] 0.219251337 0.224598930 0.229946524 0.235294118 0.240641711
## [46] 0.245989305 0.251336898 0.256684492 0.262032086 0.267379679
## [51] 0.272727273 0.278074866 0.283422460 0.288770053 0.294117647
## [56] 0.299465241 0.304812834 0.310160428 0.315508021 0.320855615
## [61] 0.326203209 0.331550802 0.336898396 0.342245989 0.347593583
## [66] 0.352941176 0.358288770 0.363636364 0.368983957 0.374331551
## [71] 0.379679144 0.385026738 0.390374332 0.395721925 0.401069519
## [76] 0.406417112 0.411764706 0.417112299 0.422459893 0.427807487
## [81] 0.433155080 0.438502674 0.443850267 0.449197861 0.454545455
## [86] 0.459893048 0.465240642 0.470588235 0.475935829 0.481283422
## [91] 0.486631016 0.491978610 0.497326203
## 
## $db
##  [1]  14.3997611  15.8955909   7.7642205   8.7050186   3.8198773
##  [6]  -1.1863576   2.6815084  -0.4667544  -1.0582274  -3.9729372
## [11]  -4.4388886  -5.3739230  -7.0574666  -3.9001118  -5.3732547
## [16]  -5.9203161  -6.0265835  -6.3073466  -8.7427682  -8.5315257
## [21]  -8.1694048  -8.8726978  -6.2228751  -5.6255282  -5.7869081
## [26]  -2.3190255  -0.6895940 -18.2510780 -12.7252753 -17.9895004
## [31] -13.1495229 -17.1809062 -12.9879824 -18.5885318 -14.0950716
## [36] -19.2823521 -16.1802599 -22.2193540 -17.1778910 -15.9864725
## [41] -13.8720977 -14.3606703 -12.3482825 -19.0159734 -12.8299009
## [46] -16.1634117 -19.2402721 -16.9940391 -19.9193164 -24.4462952
## [51] -12.2184745 -14.1347405 -11.1407485 -19.5608640 -22.3406178
## [56] -26.2948840 -15.9219942 -19.7896219 -16.3429699 -16.6545137
## [61] -13.2991307 -20.9465604 -12.5392646 -23.3015328 -21.5282804
## [66] -17.4748203 -26.6390880 -24.5923917 -17.5647166 -19.7282494
## [71] -17.6845420 -23.0825883 -19.0743556 -19.4958696 -19.5138333
## [76] -20.4270518 -18.7343253 -16.8982195 -14.1306252 -14.2054286
## [81] -15.0466299 -13.4474119 -14.7662082 -15.2848488 -14.7767269
## [86] -17.2543647 -17.5524024 -16.7680943 -23.6931393 -19.6888637
## [91] -22.9328710 -17.1805763 -15.7422487
## 
## $dbz
##  [1]  12.1744244  11.8068871  11.1912920  10.3234897   9.1987290
##  [6]   7.8134517   6.1690554   4.2794940   2.1855400  -0.0227779
## [11]  -2.1856123  -4.0930116  -5.5814407  -6.6274820  -7.3078903
## [16]  -7.6988921  -7.8526856  -7.8271468  -7.6919809  -7.5061868
## [21]  -7.3030888  -7.0977786  -6.9055053  -6.7545171  -6.6859448
## [26]  -6.7448940  -6.9710578  -7.3933729  -8.0283548  -8.8797888
## [31]  -9.9376266 -11.1745795 -12.5395304 -13.9486256 -15.2804168
## [36] -16.3913300 -17.1665125 -17.5812779 -17.7097183 -17.6695854
## [41] -17.5637030 -17.4562894 -17.3753921 -17.3239844 -17.2913960
## [46] -17.2628022 -17.2262162 -17.1767997 -17.1186074 -17.0640034
## [51] -17.0309217 -17.0382729 -17.1002705 -17.2209453 -17.3905128
## [56] -17.5857017 -17.7761429 -17.9367144 -18.0606665 -18.1649087
## [61] -18.2832201 -18.4521578 -18.6979308 -19.0286543 -19.4319317
## [66] -19.8765383 -20.3176913 -20.7053075 -20.9929698 -21.1442791
## [71] -21.1356798 -20.9585908 -20.6228760 -20.1586258 -19.6118560
## [76] -19.0349178 -18.4769322 -17.9782467 -17.5690679 -17.2703544
## [81] -17.0951532 -17.0493949 -17.1318377 -17.3332302 -17.6350503
## [86] -18.0085946 -18.4157858 -18.8133902 -19.1613062 -19.4324915
## [91] -19.6189508 -19.7292633 -19.7788462
  1. Differencing Data: much more stationary data and have surfaced a seasonaly component = 7 seen on ACF peaks at 7, 14, 21
#b
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

  1. Seaonality Transformation: stationary data, and an ACF that reflects data other than whitenoise to be modeled.
#c
inpos.diff.seas = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

  1. Model IDing: diagnose with aic.wge to determine phi, thetas: Both AIC/BIC select ARMA(4,2)
#d
aic5.wge(inpos.diff.seas)
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  aic
##       p    q        aic
## 15    4    2   15.62368
## 9     2    2   15.74701
## 18    5    2   15.74771
## 12    3    2   15.75351
## 7     2    0   15.75551
aic5.wge(inpos.diff.seas, type = "bic")
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  bic
##       p    q        bic
## 15    4    2   15.74832
## 2     0    1   15.79201
## 7     2    0   15.80893
## 4     1    0   15.81457
## 3     0    2   15.81522
  1. White Noise Test: Reject the H_null and accept this data is not white noise
#e
acf(inpos.diff.seas)

ljung.wge(inpos.diff.seas)$pval
## Obs -0.3711871 -0.02127283 0.07854862 0.05194601 -0.1044838 0.3276797 -0.4658847 0.1565306 0.1284186 -0.1060791 -0.07693327 0.1030723 -0.06391648 0.05002847 0.004471288 -0.03640442 -0.01551996 0.04133725 -0.01110058 -0.08153895 0.06570572 -0.03709739 -0.02390062 0.04588417
## [1] 1.246447e-12
  1. Estiamte Phis, Thetas
  2. Model Building
#f
est.inpos.seas = est.arma.wge(inpos.diff.seas, p = 4, q = 2)
## 
## Coefficients of Original polynomial:  
## -1.8727 -1.4974 -0.4735 0.0283 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1+1.1293B+0.6983B^2   -0.8086+-0.8822i      0.8357       0.3681
## 1+0.7944B             -1.2588               0.7944       0.5000
## 1-0.0510B              19.6244               0.0510       0.0000
##   
## 
mean(us$positiveIncrease)
## [1] 22567.12
#g
  1. Forecast
  • 7 Day
#7 day
inpos.preds7 = fore.aruma.wge(us$positiveIncrease, phi = est.inpos.seas$phi, theta = est.inpos.seas$theta, d=1, s=7, n.ahead = 7)

  • 90 Day
#90 day
inpos.preds90 = fore.aruma.wge(us$positiveIncrease, phi = est.inpos.seas$phi, theta = est.inpos.seas$theta, d=1, s=7, n.ahead = 90)

  1. Plotting forecasts
  • 7 Day Forecast
#7 day forecast
plot(seq(1,187,1), us$positiveIncrease, type = "l", xlim = c(0,195), ylim = c(0,80000), ylab = "Increase in COVID Cases", main = "7 Day Increase in COVID Cases Forecast")
lines(seq(188, 194,1), inpos.preds7$f, type = "l", col = "red")

  • 90 Day Forecast
#90 day forecast
plot(seq(1,187,1), us$positiveIncrease, type = "l", xlim = c(0,280), ylim = c(0,80000), ylab = "Increase in COVID Cases", main = "90 Day Increase in COVID Cases Forecast")
lines(seq(188, 277,1), inpos.preds90$f, type = "l", col = "red")

Forecast Dependent Variable: MLR w/ Correlated Errors for Currently Hospitalized Patients using Positive Increase

  1. Data prep and recognizing differencing and seaonal components of Currently Hospitalized
#Selecting only those dates with reported current hospitilizations
us <- dplyr::slice(us,56:187)
invisible(us)
#Stationarize Currently Hospitalized Variable
us.diff = artrans.wge(us$hospitalizedCurrently, phi.tr = 1, lag.max = 100)

acf(us.diff)
us.diff.seas = artrans.wge(us.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

acf(us.diff.seas)
#Stationarize Increase Positive Variable
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

acf(inpos.diff)
inpos.diff.seas = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

acf(inpos.diff.seas)

  1. Fit simple regression model predicting hospitalizedCurrently with positiveIncrease.
mlr.fit = lm(us.diff.seas~inpos.diff.seas, data=cbind(data.frame(us.diff.seas), data.frame(inpos.diff.seas)))
summary(mlr.fit)
## 
## Call:
## lm(formula = us.diff.seas ~ inpos.diff.seas, data = cbind(data.frame(us.diff.seas), 
##     data.frame(inpos.diff.seas)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7149.3  -554.2    73.4   611.2  6506.7 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     -13.32614  143.23299  -0.093  0.92603   
## inpos.diff.seas   0.11733    0.04227   2.776  0.00637 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1595 on 122 degrees of freedom
## Multiple R-squared:  0.0594, Adjusted R-squared:  0.0517 
## F-statistic: 7.705 on 1 and 122 DF,  p-value: 0.006375
AIC(mlr.fit)
## [1] 2184.712
acf(mlr.fit$residuals)

plot(mlr.fit$residuals)

  1. Diagnose with aic.wge and see what the p, q would be with residuals from above model
  • Below we can see that our aic.wge function has selected an ARMA(1,2) model for modeling our cmort information.
mlr.phis= aic.wge(mlr.fit$residuals)
mlr.phis
## $type
## [1] "aic"
## 
## $value
## [1] 14.53403
## 
## $p
## [1] 5
## 
## $q
## [1] 2
## 
## $phi
## [1] -0.2465069 -0.4219100  0.3119255  0.2116968  0.2024803
## 
## $theta
## [1] -0.4657537 -0.9566014
## 
## $vara
## [1] 1803071
  1. Now forcast with arima function with phi’s from above coefficients and ARMA(1,2) model,
mlr.fit.arima = arima(us$hospitalizedCurrently, order = c(5,1,2), seasonal = list(order = c(1,0,0), period = 7), xreg = us$positiveIncrease)
## Warning in log(s2): NaNs produced
AIC(mlr.fit.arima)
## [1] 2224.09
  1. Run test to see if data is white noise: confirmed white noise, continue forward
acf(mlr.fit.arima$resid) 

ltest1 = ljung.wge(mlr.fit.arima$resid) 
## Obs -0.006649913 -0.009334226 -0.01070918 0.01827832 -0.02771128 -0.01061308 -0.01123975 0.08781402 -0.05809279 0.07500195 0.05661161 -0.08527475 -0.007553036 0.1098068 0.002427647 -0.010962 -0.04242877 -0.001526285 -0.05612271 -0.03313627 0.02930685 -0.007307643 -0.06056842 -0.03342384
ltest1$pval
## [1] 0.999208
ltest2 = ljung.wge(mlr.fit.arima$resid, K= 48)
## Obs -0.006649913 -0.009334226 -0.01070918 0.01827832 -0.02771128 -0.01061308 -0.01123975 0.08781402 -0.05809279 0.07500195 0.05661161 -0.08527475 -0.007553036 0.1098068 0.002427647 -0.010962 -0.04242877 -0.001526285 -0.05612271 -0.03313627 0.02930685 -0.007307643 -0.06056842 -0.03342384 0.08407797 -0.006875273 -0.002312484 0.06777622 -0.09591675 -0.06198868 -0.01331578 -0.03441461 -0.04091228 0.0270524 0.06615092 0.06538511 -0.09261011 -0.04939741 -0.03187124 -0.04934637 -0.01834029 -0.02031817 -0.001873916 0.02337696 -0.004828854 0.006148926 -0.002455308 0.004810674
ltest2$pval
## [1] 0.9999866
  1. Load the forecasted increase positive in a data frame
  • 7 Day
#7 Day Case Increase
regs7 = data.frame(positiveIncrease = inpos.preds7$f)
invisible(regs7)
  • 90 Day
#90 Day Case Increase
regs90 = data.frame(positiveIncrease = inpos.preds90$f)
invisible(regs90)
  1. Predictions
  • 7 Day
mlr1.preds7 = predict(mlr.fit.arima, newxreg = regs7, n.ahead = 7)
invisible(mlr1.preds7)
  • 90 Day
mlr1.preds90 = predict(mlr.fit.arima, newxreg = regs90, n.ahead =90)
invisible(mlr1.preds90)
  1. Plotted Forecasts
  • 7 Day
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,140), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "7 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,139,1), mlr1.preds7$pred, type = "l", col = "red")

  • 90 Day
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,223), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "90 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,222,1), mlr1.preds90$pred, type = "l", col = "red")

  1. Windowed ASE

Forecast Dependent Variable: MLR w/ Correlated Errors for Currently Hospitalized Patients w/ Positive Increase & Trend

  1. Fit simple regression model predicting hospitalizedCurrently with positiveIncrease and trend
#creating trend
time <- seq(1,124,1)
#fitting model
mlr.fit.t = lm(us.diff.seas~inpos.diff.seas+time, data=cbind(data.frame(us.diff.seas), data.frame(inpos.diff.seas)))
summary(mlr.fit.t)
## 
## Call:
## lm(formula = us.diff.seas ~ inpos.diff.seas + time, data = cbind(data.frame(us.diff.seas), 
##     data.frame(inpos.diff.seas)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7111.0  -524.5    73.9   617.7  6540.6 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      32.35112  289.27171   0.112  0.91114   
## inpos.diff.seas   0.11724    0.04244   2.762  0.00663 **
## time             -0.73096    4.01658  -0.182  0.85590   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1601 on 121 degrees of freedom
## Multiple R-squared:  0.05966,    Adjusted R-squared:  0.04412 
## F-statistic: 3.839 on 2 and 121 DF,  p-value: 0.02419
AIC(mlr.fit.t)
## [1] 2186.678
acf(mlr.fit.t$residuals)

plot(mlr.fit.t$residuals)

  1. Diagnose with aic.wge and see what the p, q would be with residuals from above model
    1. Below we can see that our aic.wge function has selected an ARMA(5,2) model for modeling our cmort information.
mlr.phis= aic.wge(mlr.fit.t$residuals)
mlr.phis
## $type
## [1] "aic"
## 
## $value
## [1] 14.53329
## 
## $p
## [1] 5
## 
## $q
## [1] 2
## 
## $phi
## [1] -0.2466126 -0.4220817  0.3118538  0.2119376  0.2027949
## 
## $theta
## [1] -0.4655288 -0.9564592
## 
## $vara
## [1] 1801724
  1. Now forcast with arima function with phi’s from above coefficients and ARMA(1,2) model,
Time <- seq(1,132,1)
mlr.fit.arima.t = arima(us$hospitalizedCurrently, order = c(5,1,2), seasonal = list(order = c(1,0,0), period = 7), xreg = cbind(us$positiveIncrease, Time))
AIC(mlr.fit.arima.t)
## [1] 2230.756
  1. Run test to see if data is white noise; confirmed white noise continue forward.
acf(mlr.fit.arima.t$resid) 

ltest1 = ljung.wge(mlr.fit.arima.t$resid) 
## Obs 0.002267493 0.004743329 0.004085712 0.008584111 0.003613306 0.03592544 -0.02552583 0.08439651 -0.1060931 0.01276058 0.02640347 -0.06061388 0.05268146 0.1458814 0.008052735 -0.04255608 -0.0935536 -0.04208266 -0.0498013 -0.001717023 0.05934226 0.01341049 -0.07137289 -0.07279666
ltest1$pval
## [1] 0.9823161
ltest2 = ljung.wge(mlr.fit.arima.t$resid, K= 48)
## Obs 0.002267493 0.004743329 0.004085712 0.008584111 0.003613306 0.03592544 -0.02552583 0.08439651 -0.1060931 0.01276058 0.02640347 -0.06061388 0.05268146 0.1458814 0.008052735 -0.04255608 -0.0935536 -0.04208266 -0.0498013 -0.001717023 0.05934226 0.01341049 -0.07137289 -0.07279666 0.05320601 -0.01313386 0.01794322 0.1027134 -0.07595104 -0.06199621 -0.03306343 -0.06206906 -0.03944973 0.03988941 0.07572305 0.07297186 -0.1029145 -0.06857141 -0.04403521 -0.04289784 0.005428483 -0.003871671 -7.833593e-05 0.005208028 -0.03085214 -0.01036108 0.00645461 0.02868683
ltest2$pval
## [1] 0.9990337
  1. Load the forecasted increase positive in a data frame
    1. 7 Day
#7 Day Case Increase
regs7t = data.frame(cbind(positiveIncrease = inpos.preds7$f, Time = seq(133,139)))
invisible(regs7)
b. 90 Day
#90 Day Case Increase
regs90t = data.frame(cbind(positiveIncrease = inpos.preds90$f, Time = seq(133,222)))
invisible(regs90)
  1. Predictions
    1. 7 Day
mlr1.preds7.t = predict(mlr.fit.arima.t, newxreg = regs7t, n.ahead = 7)
invisible(mlr1.preds7.t)
b. 90 Day
mlr1.preds90.t = predict(mlr.fit.arima.t, newxreg = regs90t, n.ahead =90)
invisible(mlr1.preds90.t)
  1. Plotted Forecasts
    1. 7 Day
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,140), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "7 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,139,1), mlr1.preds7.t$pred, type = "l", col = "red")

b. 90 Day

plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,223), ylim = c(0,85000), ylab = "Currently Hospitalized COVID Cases", main = "90 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,222,1), mlr1.preds90.t$pred, type = "l", col = "red")

  1. Windowed ASE
  1. 7 Day: 369,785,216
phis = mlr.phis$phi
thetas = mlr.phis$theta

trainingSize = 20
horizon = 7
ASEHolder = numeric()

for( i in 1:(124-(trainingSize + horizon) + 1))
{
  forecasts = fore.aruma.wge(us$hospitalizedCurrently[i:(i+(trainingSize-1))],phi = phis, theta = thetas, s = 7, d = 1,n.ahead = horizon)
  
  ASE = mean((us$hospitalizedCurrently[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - mlr1.preds7.t$pred)^2)
         
  ASEHolder[i] = ASE
}
ASEHolder
##  [1] 156784097  81468394  46794011  20158053  10075549   4247909   2536964
##  [8]   1954887   2724789   3478851   3398630   2537931   2240228   2141373
## [15]   2922274   4328430   5618462   6902646   8755834  10550858  13373208
## [22]  17834238  22852046  26712915  31252713  36683951  45986421  57837302
## [29]  73104265  90379284 105247611 120809992 137709830 155968235 175820326
## [36] 198604255 219873583 243790642 265573514 285773904 307106009 331522863
## [43] 354469460 377747700 400740703 419259524 436308656 453354771 471097774
## [50] 489807847 512541593 538652070 572796462 608070438 639639595 667002781
## [57] 692624259 715931788 736840013 755442379 780285643 806871492 835829489
## [64] 865656328 894120241 918750075 937836561 948640485 952994160 957055984
## [71] 954153568 946215541 929638610 908358035 879601578 853181389 815665729
## [78] 780725150 742790976 705231840 666766187 626782470 584254960 545301500
## [85] 503155266 463843112 425716302 388787951 354469722 297543457 241934983
## [92] 187058113 136636337  95276468  59832810  28643460  20500919  13079159
hist(ASEHolder)

WindowedASE = mean(ASEHolder)
summary(ASEHolder)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##   1954887  32610523 302324733 369785216 666943633 957055984
WindowedASE
## [1] 369785216
  1. 90 Day: 1,054,226,906
phis = mlr.phis$phi
thetas = mlr.phis$theta

trainingSize = 20
horizon = 90
ASEHolder = numeric()

for( i in 1:(124-(trainingSize + horizon) + 1))
{
  forecasts = fore.aruma.wge(us$hospitalizedCurrently[i:(i+(trainingSize-1))],phi = phis, theta = thetas, s = 7, d = 1,n.ahead = horizon)
  
  ASE = mean((us$hospitalizedCurrently[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - mlr1.preds90.t$pred)^2)
         
  ASEHolder[i] = ASE
}
ASEHolder
##  [1] 1074315371 1074565667 1076796973 1077260763 1077296084 1076804189
##  [7] 1070559599 1064378347 1057995663 1050717725 1042087440 1032810213
## [13] 1022673853 1012502262 1002639450
hist(ASEHolder)

WindowedASE = mean(ASEHolder)
summary(ASEHolder)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.003e+09 1.037e+09 1.064e+09 1.054e+09 1.076e+09 1.077e+09
WindowedASE
## [1] 1054226906

Forecast Dependent Variable: MLR w/ Correlated Errors (lagged variables) for Currently Hospitalized Patients w/ Positive Increase & Trend

With a quick check, we can see that there is no lag correlationed between the increase of COVID patients and hospitalized patients, like we theorized there might be. So we will not model an MLR w/ correlated errors on lagged variables.

ccf(us$positiveIncrease,us$hospitalizedCurrently)

Multivariate VAR Model

Since we are working with a seasonal and transformation component but it requires an additional transformation for the total positive COVID cases to become stationary for evaluation, we will only use positive increase of cases to predict currently hospitzliaed cases for the VAR model.

  1. Differenced and Seasonal Transformation VAR Model
#Positive Cases Transformations
#pos.diff = artrans.wge(us$positive, phi.tr = 1, lag.max = 100)
#pos.diff.2 = artrans.wge(pos.diff, phi.tr = 1, lag.max = 100)
#pos.trans = artrans.wge(pos.diff.2,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)
#Increase Positive Transformation
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

inpos.trans = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

#Current Hospitalization Transformations
us.diff = artrans.wge(us$hospitalizedCurrently, phi.tr = 1, lag.max = 100)

currhosp.trans = artrans.wge(us.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

  1. Use Var select to find best model and fit model: p = 8 for lowest AIC
#VARselect to choose lag
VARselect(cbind(currhosp.trans, inpos.trans), lag.max = 10, type= "both")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      7      7      2      7 
## 
## $criteria
##                   1            2            3            4            5
## AIC(n) 3.093857e+01 3.083329e+01 3.079091e+01 3.076976e+01 3.069554e+01
## HQ(n)  3.101650e+01 3.095018e+01 3.094676e+01 3.096458e+01 3.092932e+01
## SC(n)  3.113058e+01 3.112131e+01 3.117493e+01 3.124979e+01 3.127158e+01
## FPE(n) 2.731960e+13 2.459304e+13 2.357880e+13 2.309552e+13 2.145770e+13
##                   6            7            8            9           10
## AIC(n) 3.066656e+01 3.051427e+01 3.053424e+01 3.059674e+01 3.059649e+01
## HQ(n)  3.093930e+01 3.082598e+01 3.088492e+01 3.098638e+01 3.102509e+01
## SC(n)  3.133860e+01 3.128233e+01 3.139830e+01 3.155681e+01 3.165257e+01
## FPE(n) 2.086403e+13 1.793906e+13 1.833019e+13 1.955149e+13 1.959499e+13
#fit model
usdiff.var.fit = VAR(cbind(currhosp.trans, inpos.trans), type = "both", p = 2)
  1. Predictions for Difference
    1. 7 Day
#7 day predictions
pred.var7 = predict(usdiff.var.fit, n.ahead = 7)
pred.var7$fcst$currhosp.trans
##            fcst     lower    upper       CI
## [1,]  -42.12350 -2942.410 2858.163 2900.286
## [2,] -385.17371 -3346.493 2576.145 2961.319
## [3,]  -80.65766 -3262.529 3101.214 3181.871
## [4,] -124.14161 -3327.324 3079.041 3203.182
## [5,]  -46.62794 -3284.538 3191.282 3237.910
## [6,]  -43.58380 -3288.207 3201.039 3244.623
## [7,]  -11.23201 -3262.277 3239.813 3251.045
b. 90 Day
#90 day predictions
pred.var90 = predict(usdiff.var.fit, n.ahead = 90)
pred.var90$fcst$currhosp.trans
##              fcst     lower    upper       CI
##  [1,]  -42.123501 -2942.410 2858.163 2900.286
##  [2,] -385.173707 -3346.493 2576.145 2961.319
##  [3,]  -80.657665 -3262.529 3101.214 3181.871
##  [4,] -124.141611 -3327.324 3079.041 3203.182
##  [5,]  -46.627942 -3284.538 3191.282 3237.910
##  [6,]  -43.583795 -3288.207 3201.039 3244.623
##  [7,]  -11.232008 -3262.277 3239.813 3251.045
##  [8,]   -6.349200 -3259.226 3246.528 3252.877
##  [9,]    7.305513 -3246.870 3261.481 3254.175
## [10,]   12.421029 -3242.210 3267.053 3254.632
## [11,]   18.773716 -3236.131 3273.679 3254.905
## [12,]   22.480950 -3232.533 3277.495 3255.014
## [13,]   26.203504 -3228.870 3281.277 3255.073
## [14,]   28.923414 -3226.175 3284.022 3255.099
## [15,]   31.503265 -3223.608 3286.615 3255.112
## [16,]   33.699320 -3221.418 3288.817 3255.117
## [17,]   35.776726 -3219.344 3290.897 3255.120
## [18,]   37.697587 -3217.424 3292.819 3255.122
## [19,]   39.549962 -3215.572 3294.672 3255.122
## [20,]   41.334045 -3213.789 3296.457 3255.123
## [21,]   43.082251 -3212.041 3298.205 3255.123
## [22,]   44.799892 -3210.323 3299.923 3255.123
## [23,]   46.499486 -3208.623 3301.622 3255.123
## [24,]   48.185069 -3206.938 3303.308 3255.123
## [25,]   49.861824 -3205.261 3304.985 3255.123
## [26,]   51.532055 -3203.591 3306.655 3255.123
## [27,]   53.198022 -3201.925 3308.321 3255.123
## [28,]   54.860928 -3200.262 3309.984 3255.123
## [29,]   56.521787 -3198.601 3311.645 3255.123
## [30,]   58.181203 -3196.942 3313.304 3255.123
## [31,]   59.839642 -3195.283 3314.963 3255.123
## [32,]   61.497398 -3193.626 3316.620 3255.123
## [33,]   63.154688 -3191.968 3318.278 3255.123
## [34,]   64.811654 -3190.311 3319.935 3255.123
## [35,]   66.468399 -3188.655 3321.591 3255.123
## [36,]   68.124990 -3186.998 3323.248 3255.123
## [37,]   69.781476 -3185.341 3324.904 3255.123
## [38,]   71.437889 -3183.685 3326.561 3255.123
## [39,]   73.094252 -3182.029 3328.217 3255.123
## [40,]   74.750580 -3180.372 3329.874 3255.123
## [41,]   76.406884 -3178.716 3331.530 3255.123
## [42,]   78.063172 -3177.060 3333.186 3255.123
## [43,]   79.719449 -3175.404 3334.842 3255.123
## [44,]   81.375717 -3173.747 3336.499 3255.123
## [45,]   83.031981 -3172.091 3338.155 3255.123
## [46,]   84.688240 -3170.435 3339.811 3255.123
## [47,]   86.344498 -3168.778 3341.467 3255.123
## [48,]   88.000753 -3167.122 3343.124 3255.123
## [49,]   89.657007 -3165.466 3344.780 3255.123
## [50,]   91.313260 -3163.810 3346.436 3255.123
## [51,]   92.969513 -3162.153 3348.092 3255.123
## [52,]   94.625766 -3160.497 3349.749 3255.123
## [53,]   96.282018 -3158.841 3351.405 3255.123
## [54,]   97.938270 -3157.185 3353.061 3255.123
## [55,]   99.594521 -3155.528 3354.717 3255.123
## [56,]  101.250773 -3153.872 3356.374 3255.123
## [57,]  102.907025 -3152.216 3358.030 3255.123
## [58,]  104.563276 -3150.560 3359.686 3255.123
## [59,]  106.219528 -3148.903 3361.342 3255.123
## [60,]  107.875779 -3147.247 3362.999 3255.123
## [61,]  109.532031 -3145.591 3364.655 3255.123
## [62,]  111.188282 -3143.935 3366.311 3255.123
## [63,]  112.844534 -3142.278 3367.967 3255.123
## [64,]  114.500785 -3140.622 3369.624 3255.123
## [65,]  116.157037 -3138.966 3371.280 3255.123
## [66,]  117.813288 -3137.310 3372.936 3255.123
## [67,]  119.469540 -3135.653 3374.592 3255.123
## [68,]  121.125791 -3133.997 3376.249 3255.123
## [69,]  122.782043 -3132.341 3377.905 3255.123
## [70,]  124.438294 -3130.685 3379.561 3255.123
## [71,]  126.094546 -3129.028 3381.218 3255.123
## [72,]  127.750797 -3127.372 3382.874 3255.123
## [73,]  129.407049 -3125.716 3384.530 3255.123
## [74,]  131.063300 -3124.060 3386.186 3255.123
## [75,]  132.719552 -3122.403 3387.843 3255.123
## [76,]  134.375803 -3120.747 3389.499 3255.123
## [77,]  136.032055 -3119.091 3391.155 3255.123
## [78,]  137.688306 -3117.435 3392.811 3255.123
## [79,]  139.344558 -3115.778 3394.468 3255.123
## [80,]  141.000809 -3114.122 3396.124 3255.123
## [81,]  142.657061 -3112.466 3397.780 3255.123
## [82,]  144.313312 -3110.810 3399.436 3255.123
## [83,]  145.969564 -3109.153 3401.093 3255.123
## [84,]  147.625815 -3107.497 3402.749 3255.123
## [85,]  149.282067 -3105.841 3404.405 3255.123
## [86,]  150.938318 -3104.185 3406.061 3255.123
## [87,]  152.594570 -3102.528 3407.718 3255.123
## [88,]  154.250821 -3100.872 3409.374 3255.123
## [89,]  155.907073 -3099.216 3411.030 3255.123
## [90,]  157.563324 -3097.560 3412.686 3255.123
  1. Calculate Actual Forecasts from Predicted Differences
    1. 7 Day
startingPoints7 = us$hospitalizedCurrently[126:132]
currHospForecasts7 = pred.var7$fcst$currhosp.trans[,1:3] + startingPoints7
b. 90 Day
startingPoints90 = us$hospitalizedCurrently[43:132]
currHospForecasts90 = pred.var90$fcst$currhosp.trans[,1:3] + startingPoints90
  1. Plotting Forecasts
    1. 7 Day
#7 day Forecasts
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l", xlim = c(0,139), ylab = "Currently Hospitalized COVID Patients", main = "7 Day Currently Hospitalized Patients Forecast")
lines(seq(133,139,1), as.data.frame(currHospForecasts7)$fcst, type = "l", col = "red")

b. 90 Day

#90 day Forecasts
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l", xlim = c(0,222), ylab = "Currently Hospitalized COVID Patients", main = "90 Day Currently Hospitalized Patients Forecast")
lines(seq(133,222,1), as.data.frame(currHospForecasts90)$fcst, type = "l", col = "red")

  1. Calculate ASE # WHERE YOU LEFT OFF
#varASE = mean((us$hospitalizedCurrently[123:132]-pred.var$fcst$y1[1:10])^2)
#varASE
  1. Windowed ASE

Multilayer Perceptron Model (MLP) / NN

Since we have been looking at additional regressors for all our above models and know that the best regressor to leverage will be positive increase of COVID cases. This regressor reflects similar behavior to that of current hospitlizations and can be model properly against hospitalizations. 1. Create Current Hospitalized ts variable

tUScurrhop = ts(us$hospitalizedCurrently)
  1. Create data frame of regressor: positive increase covid cases variable
tUSx = data.frame(positiveIncrease = ts(us$positiveIncrease))
  1. Forecast of positive increase of COVID cases
    1. 7 Day Forecast for Regressor
fit.mlp.new = mlp(ts(tUSx), reps = 50, comb = "mean")
mlp.new.fore7 = forecast(fit.mlp.new, h = 7)
mlp.new.fore7
##     Point Forecast
## 133       60022.23
## 134       63703.79
## 135       65465.01
## 136       66514.70
## 137       65662.89
## 138       65030.57
## 139       64549.13
b. 90 Day Forecast for Regressor
mlp.new.fore90 = forecast(fit.mlp.new, h = 90)
mlp.new.fore90
##     Point Forecast
## 133       60022.23
## 134       63703.79
## 135       65465.01
## 136       66514.70
## 137       65662.89
## 138       65030.57
## 139       64549.13
## 140       64841.09
## 141       65153.74
## 142       65416.89
## 143       65407.32
## 144       65312.68
## 145       65086.79
## 146       65086.33
## 147       65141.91
## 148       65285.80
## 149       65337.72
## 150       65339.22
## 151       65203.62
## 152       65132.75
## 153       65136.06
## 154       65237.91
## 155       65310.28
## 156       65340.96
## 157       65272.21
## 158       65162.37
## 159       65138.63
## 160       65197.44
## 161       65281.46
## 162       65327.43
## 163       65312.92
## 164       65202.48
## 165       65154.10
## 166       65163.65
## 167       65246.95
## 168       65304.72
## 169       65326.67
## 170       65256.25
## 171       65172.70
## 172       65153.24
## 173       65212.76
## 174       65280.81
## 175       65320.37
## 176       65297.83
## 177       65197.72
## 178       65161.98
## 179       65180.20
## 180       65253.64
## 181       65302.96
## 182       65317.58
## 183       65239.50
## 184       65176.04
## 185       65161.42
## 186       65224.66
## 187       65283.02
## 188       65317.72
## 189       65283.32
## 190       65192.38
## 191       65163.94
## 192       65194.32
## 193       65261.09
## 194       65305.04
## 195       65309.56
## 196       65223.43
## 197       65175.53
## 198       65169.22
## 199       65235.08
## 200       65287.03
## 201       65316.29
## 202       65267.72
## 203       65188.04
## 204       65164.37
## 205       65207.13
## 206       65268.25
## 207       65308.30
## 208       65300.29
## 209       65210.14
## 210       65173.17
## 211       65178.87
## 212       65244.87
## 213       65292.20
## 214       65314.14
## 215       65250.46
## 216       65184.51
## 217       65165.45
## 218       65218.85
## 219       65274.88
## 220       65311.36
## 221       65288.90
## 222       65200.09
  1. Combine observed new cases + forecast new cases
    1. 7 Day regressor var
new.regressor7 <- data.frame(c(us$positiveIncrease, mlp.new.fore7$mean))
new.regressor7 
##     c.us.positiveIncrease..mlp.new.fore7.mean.
## 1                                      3613.00
## 2                                      3171.00
## 3                                      4671.00
## 4                                      6255.00
## 5                                      6885.00
## 6                                      9259.00
## 7                                     11442.00
## 8                                     10632.00
## 9                                     12873.00
## 10                                    17656.00
## 11                                    19044.00
## 12                                    19724.00
## 13                                    19655.00
## 14                                    21897.00
## 15                                    24694.00
## 16                                    25773.00
## 17                                    27992.00
## 18                                    31945.00
## 19                                    33252.00
## 20                                    25550.00
## 21                                    28937.00
## 22                                    30737.00
## 23                                    30534.00
## 24                                    34419.00
## 25                                    34351.00
## 26                                    30561.00
## 27                                    27844.00
## 28                                    25133.00
## 29                                    25631.00
## 30                                    30298.00
## 31                                    30923.00
## 32                                    31964.00
## 33                                    27963.00
## 34                                    27468.00
## 35                                    25872.00
## 36                                    26333.00
## 37                                    28916.00
## 38                                    31784.00
## 39                                    34174.00
## 40                                    35901.00
## 41                                    27380.00
## 42                                    22605.00
## 43                                    25219.00
## 44                                    26641.00
## 45                                    29568.00
## 46                                    33056.00
## 47                                    29185.00
## 48                                    25767.00
## 49                                    22523.00
## 50                                    22449.00
## 51                                    25056.00
## 52                                    27490.00
## 53                                    27605.00
## 54                                    24810.00
## 55                                    21504.00
## 56                                    18236.00
## 57                                    22663.00
## 58                                    21193.00
## 59                                    26705.00
## 60                                    24710.00
## 61                                    24701.00
## 62                                    20171.00
## 63                                    20992.00
## 64                                    20853.00
## 65                                    21319.00
## 66                                    26513.00
## 67                                    24555.00
## 68                                    21590.00
## 69                                    20105.00
## 70                                    18798.00
## 71                                    16629.00
## 72                                    19385.00
## 73                                    22601.00
## 74                                    23522.00
## 75                                    23762.00
## 76                                    21639.00
## 77                                    20415.00
## 78                                    19982.00
## 79                                    20312.00
## 80                                    20839.00
## 81                                    23352.00
## 82                                    23106.00
## 83                                    18823.00
## 84                                    17012.00
## 85                                    17175.00
## 86                                    20803.00
## 87                                    22032.00
## 88                                    23477.00
## 89                                    25341.00
## 90                                    21373.00
## 91                                    18510.00
## 92                                    23435.00
## 93                                    23873.00
## 94                                    27527.00
## 95                                    31046.00
## 96                                    31994.00
## 97                                    27284.00
## 98                                    27017.00
## 99                                    33021.00
## 100                                   38684.00
## 101                                   39072.00
## 102                                   44421.00
## 103                                   43783.00
## 104                                   41857.00
## 105                                   39175.00
## 106                                   47462.00
## 107                                   50674.00
## 108                                   53826.00
## 109                                   54223.00
## 110                                   54734.00
## 111                                   45789.00
## 112                                   41600.00
## 113                                   51766.00
## 114                                   62147.00
## 115                                   58836.00
## 116                                   66645.00
## 117                                   63007.00
## 118                                   60978.00
## 119                                   58465.00
## 120                                   62879.00
## 121                                   65382.00
## 122                                   70953.00
## 123                                   77233.00
## 124                                   65180.00
## 125                                   64884.00
## 126                                   56971.00
## 127                                   63642.00
## 128                                   69150.00
## 129                                   71027.00
## 130                                   75193.00
## 131                                   65413.00
## 132                                   61713.00
## 133                                   60022.23
## 134                                   63703.79
## 135                                   65465.01
## 136                                   66514.70
## 137                                   65662.89
## 138                                   65030.57
## 139                                   64549.13
b. 90 day regressor var
new.regressor90 <- data.frame(c(us$positiveIncrease, mlp.new.fore90$mean))
new.regressor90
##     c.us.positiveIncrease..mlp.new.fore90.mean.
## 1                                       3613.00
## 2                                       3171.00
## 3                                       4671.00
## 4                                       6255.00
## 5                                       6885.00
## 6                                       9259.00
## 7                                      11442.00
## 8                                      10632.00
## 9                                      12873.00
## 10                                     17656.00
## 11                                     19044.00
## 12                                     19724.00
## 13                                     19655.00
## 14                                     21897.00
## 15                                     24694.00
## 16                                     25773.00
## 17                                     27992.00
## 18                                     31945.00
## 19                                     33252.00
## 20                                     25550.00
## 21                                     28937.00
## 22                                     30737.00
## 23                                     30534.00
## 24                                     34419.00
## 25                                     34351.00
## 26                                     30561.00
## 27                                     27844.00
## 28                                     25133.00
## 29                                     25631.00
## 30                                     30298.00
## 31                                     30923.00
## 32                                     31964.00
## 33                                     27963.00
## 34                                     27468.00
## 35                                     25872.00
## 36                                     26333.00
## 37                                     28916.00
## 38                                     31784.00
## 39                                     34174.00
## 40                                     35901.00
## 41                                     27380.00
## 42                                     22605.00
## 43                                     25219.00
## 44                                     26641.00
## 45                                     29568.00
## 46                                     33056.00
## 47                                     29185.00
## 48                                     25767.00
## 49                                     22523.00
## 50                                     22449.00
## 51                                     25056.00
## 52                                     27490.00
## 53                                     27605.00
## 54                                     24810.00
## 55                                     21504.00
## 56                                     18236.00
## 57                                     22663.00
## 58                                     21193.00
## 59                                     26705.00
## 60                                     24710.00
## 61                                     24701.00
## 62                                     20171.00
## 63                                     20992.00
## 64                                     20853.00
## 65                                     21319.00
## 66                                     26513.00
## 67                                     24555.00
## 68                                     21590.00
## 69                                     20105.00
## 70                                     18798.00
## 71                                     16629.00
## 72                                     19385.00
## 73                                     22601.00
## 74                                     23522.00
## 75                                     23762.00
## 76                                     21639.00
## 77                                     20415.00
## 78                                     19982.00
## 79                                     20312.00
## 80                                     20839.00
## 81                                     23352.00
## 82                                     23106.00
## 83                                     18823.00
## 84                                     17012.00
## 85                                     17175.00
## 86                                     20803.00
## 87                                     22032.00
## 88                                     23477.00
## 89                                     25341.00
## 90                                     21373.00
## 91                                     18510.00
## 92                                     23435.00
## 93                                     23873.00
## 94                                     27527.00
## 95                                     31046.00
## 96                                     31994.00
## 97                                     27284.00
## 98                                     27017.00
## 99                                     33021.00
## 100                                    38684.00
## 101                                    39072.00
## 102                                    44421.00
## 103                                    43783.00
## 104                                    41857.00
## 105                                    39175.00
## 106                                    47462.00
## 107                                    50674.00
## 108                                    53826.00
## 109                                    54223.00
## 110                                    54734.00
## 111                                    45789.00
## 112                                    41600.00
## 113                                    51766.00
## 114                                    62147.00
## 115                                    58836.00
## 116                                    66645.00
## 117                                    63007.00
## 118                                    60978.00
## 119                                    58465.00
## 120                                    62879.00
## 121                                    65382.00
## 122                                    70953.00
## 123                                    77233.00
## 124                                    65180.00
## 125                                    64884.00
## 126                                    56971.00
## 127                                    63642.00
## 128                                    69150.00
## 129                                    71027.00
## 130                                    75193.00
## 131                                    65413.00
## 132                                    61713.00
## 133                                    60022.23
## 134                                    63703.79
## 135                                    65465.01
## 136                                    66514.70
## 137                                    65662.89
## 138                                    65030.57
## 139                                    64549.13
## 140                                    64841.09
## 141                                    65153.74
## 142                                    65416.89
## 143                                    65407.32
## 144                                    65312.68
## 145                                    65086.79
## 146                                    65086.33
## 147                                    65141.91
## 148                                    65285.80
## 149                                    65337.72
## 150                                    65339.22
## 151                                    65203.62
## 152                                    65132.75
## 153                                    65136.06
## 154                                    65237.91
## 155                                    65310.28
## 156                                    65340.96
## 157                                    65272.21
## 158                                    65162.37
## 159                                    65138.63
## 160                                    65197.44
## 161                                    65281.46
## 162                                    65327.43
## 163                                    65312.92
## 164                                    65202.48
## 165                                    65154.10
## 166                                    65163.65
## 167                                    65246.95
## 168                                    65304.72
## 169                                    65326.67
## 170                                    65256.25
## 171                                    65172.70
## 172                                    65153.24
## 173                                    65212.76
## 174                                    65280.81
## 175                                    65320.37
## 176                                    65297.83
## 177                                    65197.72
## 178                                    65161.98
## 179                                    65180.20
## 180                                    65253.64
## 181                                    65302.96
## 182                                    65317.58
## 183                                    65239.50
## 184                                    65176.04
## 185                                    65161.42
## 186                                    65224.66
## 187                                    65283.02
## 188                                    65317.72
## 189                                    65283.32
## 190                                    65192.38
## 191                                    65163.94
## 192                                    65194.32
## 193                                    65261.09
## 194                                    65305.04
## 195                                    65309.56
## 196                                    65223.43
## 197                                    65175.53
## 198                                    65169.22
## 199                                    65235.08
## 200                                    65287.03
## 201                                    65316.29
## 202                                    65267.72
## 203                                    65188.04
## 204                                    65164.37
## 205                                    65207.13
## 206                                    65268.25
## 207                                    65308.30
## 208                                    65300.29
## 209                                    65210.14
## 210                                    65173.17
## 211                                    65178.87
## 212                                    65244.87
## 213                                    65292.20
## 214                                    65314.14
## 215                                    65250.46
## 216                                    65184.51
## 217                                    65165.45
## 218                                    65218.85
## 219                                    65274.88
## 220                                    65311.36
## 221                                    65288.90
## 222                                    65200.09
  1. Fit model for currently hospitalized w/ regressor of new cases
mlp.fit1 = mlp(tUScurrhop, xreg = tUSx, outplot = T, comb = "mean")

plot(mlp.fit1)

  1. Forecast w/ known Positive Increase Case Variable
    1. Currently Hospitalized 7 Day Forecast w/ Positive Increaser Regressor
currhosp.fore7 = forecast(mlp.fit1, h = 7, xreg = new.regressor7)
plot(currhosp.fore7)

b. Currently Hospitalized 90 Day Forecast w/ Positive Increaser Regressor

currhosp.fore90 = forecast(mlp.fit1, h = 90, xreg = new.regressor90)
plot(currhosp.fore90)

  1. Simple ASE
ASEmlr = mean((us$hospitalizedCurrently[126:132]-currhosp.fore7$mean)^2)
ASEmlr
## [1] 12199114

MLP NN Model Analysis

We completed a default neural network model. With so many opportunities for how to actually tune neural network model we knew this would not be our best model in this case. So we moved forward with a hypertuned neural network model that allows us to calculate many windowed ASEs and compare those model against eachother.

Ensemble Model / Hypertuned NN Model

We have leveraged the tswgewrapper code above to produce a hyperparameter tuned NN model for our ensemble model. This model will work as our higher functioning ensemble model. Below are the steps to surface the best hyperparameters based on our grid search

#if (!require(devtools)) {
#  install.packages("devtools")
#}
#devtools::install_github("josephsdavid/tswgewrapped", build_vignettes = TRUE)
library(tswgewrapped)
  1. Train/ Test Data Sets
data_train.m <- data.frame(hospitalizedCurrently = us$hospitalizedCurrently[1:122], positiveIncrease = us$positiveIncrease[1:122])
data_test.m <- data.frame(hospitalizedCurrently = us$hospitalizedCurrently[123:132], positiveIncrease = us$positiveIncrease[123:132])
  1. Hyper tune parameters
# search for best NN hyperparameters in given grid
model.m = tswgewrapped::ModelBuildNNforCaret$new(data = data_train.m, var_interest = "hospitalizedCurrently",
                                               search = 'random', tuneLength = 5, parallel = TRUE,
                                               batch_size = 50, h = 7, m = 7,
                                               verbose = 1)
## Registered S3 method overwritten by 'lava':
##   method     from   
##   print.pcor greybox
## Aggregating results
## Selecting tuning parameters
## Fitting reps = 10, hd = 2, allow.det.season = TRUE on full training set
## - Total Time for training: : 52.316 sec elapsed
  1. The windowed ASEs associated with the grid of hyperparameters is shown in the table and heatmap below.
res.m <- model.m$summarize_hyperparam_results()
res.m
##   reps hd allow.det.season     RMSE      ASE   RMSESD    ASESD
## 1   10  2             TRUE 3103.429 13970241 2184.689 16296534
## 2   11  3            FALSE 3129.341 14942705 2380.109 19673017
## 3   16  1            FALSE 3559.638 16480775 2047.127 15988977
## 4   16  3             TRUE 3204.752 15391189 2373.358 19509361
## 5   22  3            FALSE 3268.191 16196849 2463.200 19979406
model.m$plot_hyperparam_results()

  1. Best Parameters shown in below table. The best hyperparameters based on this grid search are 10 repetitions and 2 hidden layers, and allow.det.season = TRUE .
best.m <- model.m$summarize_best_hyperparams()
best.m
##   reps hd allow.det.season
## 1   10  2             TRUE
  1. Windowed ASE of 13,970,241.
final.ase.m <- dplyr::filter(res.m, reps == best.m$reps &
                    hd == best.m$hd &
                    allow.det.season == best.m$allow.det.season)[['ASE']]
final.ase.m
## [1] 13970241
  1. Ensemble model characteristics and plot
# Final Model
caret_model.m = model.m$get_final_models(subset = 'a')
caret_model.m$finalModel
## MLP fit with 2 hidden nodes and 10 repetitions.
## Series modelled in differences: D1.
## Univariate lags: (1,2,3)
## 1 regressor included.
## - Regressor 1 lags: (1,3)
## Forecast combined using the median operator.
## MSE: 1065689.5459.
# Plot Final Model
plot(caret_model.m$finalModel)

  1. Forecasts
  • Model final with best hyper parameters from above and regressor of positive increaser of COVID cases.
#Ensemble model
ensemble.mlp = mlp(tUScurrhop, xreg = tUSx, outplot = T, reps = 10, hd = 2, allow.det.season = T)

ensemble.mlp
## MLP fit with 2 hidden nodes and 10 repetitions.
## Series modelled in differences: D1.
## Univariate lags: (1,2,3)
## 1 regressor included.
## - Regressor 1 lags: (1,3)
## Forecast combined using the median operator.
## MSE: 1025155.9706.
#Plot ensemble model
plot(ensemble.mlp)

  • 7 Day
fore7.m = forecast(ensemble.mlp, xreg = new.regressor7, h=7)
plot(fore7.m)

  • 90 Day
fore90.m = forecast(ensemble.mlp, xreg = new.regressor90, h=90)
plot(fore90.m)